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1. Introduction 



The recent reports of Bose-Einstein condensation 
(BEC) in weakly interacting trapped alkali gases [1-3] 
has confirmed a property of bosons first predicted in 
1924 by Bose [4] for photons, and in 1925 by Einstein 
[5] for atoms. The production of such condensates has 
opened the possibility of a new generation of atomic 
physics experiments on meso- or macro-scopic assem- 
blies of atoms in the same quantum state. 



The pure BEC phenomenon is manifested in the 
quantum statistical mechanics of noninteracting Bose 
particles; when encountered in nature it is modified by 
the effects of particle interaction. Such effects are quite 
severe in the case of superfluid "^He, which heretofore 
has been thought of as the canonical example of BEC: 
the strong interactions between "^He atoms allow only 
about 1 % of the atoms to occupy the condensate state 
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[6] . However, it has recently been shown [7] that on the 
periphery of "^He droplets, the condensate fraction ap- 
proaches unity; this is because the atomic density van- 
ishes at the extremity of the droplet, so the atom-atom 
interaction energy also tends to zero. The recent alkali 
BECs are distinguished from liquid "^He in that the alkali 
atomic density is subject to wide-ranging experimental 
control, and has been indeed taken into the regime in 
which a dilute gas approximation is applicable, and in 
which it is reasonable to think of all atoms as occupying 
the same quantum state, as in the original concept of 
BEC. However, as will be seen below, the condensate 
state is much different from that which would describe 
a noninteracting gas. It is modified by effects of interac- 
tions that are encapsulated in a mean-field description, 
which can be thought of as the Hartree approximation to 
the wavefunction of a system of Bose particles. The 
quantitative treatment of such systems is the subject of 
this paper. 

Before proceeding, we wish to note three additional 
aspects of the alkali systems that contrast with the case 
of liquid "^He, as these make the general approach to the 
problem somewhat different from those traditionally 
used to treat superfluid systems. 

First, the alkalis are confined by an external potential 
(a magnetic field or combination of magnetic and light 
fields), so their density is inhomogeneous. Thus the 
alkali BECs cannot be described adequately by a spa- 
tially uniform condensate wavefunction such as that 
which is used to describe bulk liquid "^He. Not only must 
quantitative modelling methods be modified to treat in- 
homogeneous vs homogeneous BECs, but there are 
qualitative differences as well: for negative scattering 
lengths, a small long-lived BEC can exist in the inhomo- 
geneous case [8,9], but not in a homogeneous system 
[10]. 

Second, as discussed by Cornell [1 1], the alkali BECs 
are intrinsically metastable . The equilibrium state of a 
confined alkali system at sub-microKelvin temperatures 
is a solid. However, the time for recombination of the gas 
is very long in the dilute limit, and is at least of the order 
of seconds in the current experiments. 

Third, as shown elsewhere in this Special Issue [12- 
14], the ultracold collision physics of alkali BECs is 
exceedingly complex. Although the effects of collisions 
can be encapsulated in a few parameters (scattering 
lengths), the quantitative determination of these parame- 
ters is quite difficult and remains an active area of re- 
search. The work described in this paper utilizes these 
parameters as basic input, and it should be kept in mind 
that their values are subject to significant uncertainties, 
in no case less than 10 %. 

This paper presents a partial review of work we have 
undertaken to date in the field of modeling the alkali 



BECs. This work has been based on a zero-temperature, 
mean-field formulation of the quantum mechanics of an 
externally confined system of weakly interacting Bose 
particles. Many of the results of this theory, such as 
condensate geometries, lifetimes, and excitation fre- 
quencies, can be directly compared with the data of 
current experiments. We shall use this comparison to 
assess the validity of applying mean-field theory (MET) 
to the current crop of experimental Bose-Einstein con- 
densates (BECs). 

The zero-temperature MET equations presented here 
were first derived by Bogoliubov [15] many years ago in 
order to study the superfluid "^He. The system to which 
they apply is assumed to be a weakly interacting, dilute 
gas of identical bosons, which, as noted above, does not 
provide a good description of liquid helium. However, it 
seems to fit the conditions present in a system of mag- 
netically trapped gas of neutral alkali atoms. We empha- 
size that the previous statement should not be taken to be 
true a priori^ but rather must be subjected to stringent 
experimental tests. We present here an overview of the 
comparison of MET predictions with experiment. 

The plan of the paper is as follows. In Sec. 2 we 
present a derivation of the Gross -Pi teavskii (GP) and 
Bogoliubov equations (which we have been calling here 
the "MET" equations). As part of the discussion we 
shall attempt to provide a detailed description of all of 
the approximations made in arriving at the MET equa- 
tions. In Sec. 3 we present the results of solving these 
equations for cases where comparison with experiment 
is possible. Sections 4 and 5 describe the algorithms and 
numerical procedures we have used to obtain the results 
presented in this paper and in previous work cited 
therein. Actual solution of the MET equations for cases 
of specific experimental interest is a subject that has 
developed quite recently, and we believe there is consid- 
erable scope for enhancement of computational effi- 
ciency over that attained in current practice. Such en- 
hancements will certainly be needed to go beyond the 
zero-temperature MET description of BEC. Thus the 
material in Sees. 4 and 5 is presented at a level of detail 
needed to document our approach for use by those who 
can improve upon it. 



2. Mean-Field Theory: Approximations 
and Derivations 

In this section we present a somewhat detailed deriva- 
tion of the basic zero- temperature MET equations. 
These equations consist of the Gross -Pi taevskii equa- 
tion, which describes the properties of the condensed 
part of the trapped atomic cloud, and the Bogoliubov 
equations, which describe properties of the non- 
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condensed part. We shall present two derivations of the 
MFT equations. The first derivation uses a Bogoliubov 
transformation to cast the grand-canonical hamiltonian 
for a collection of interacting bosons into the form of a 
collection of noninteracting quasi-particles with the 
condensate becoming the vacuum state. The second 
derivation uses linear-response theory [16] performed 
on the time-dependent Gross-Pitaevskii equation (which 
is itself derived from a variational principle) to obtain 
the basic MFT equations. Before presenting these 
derivations, we shall first discuss the fundamental ap- 
proximations made in modeling a cloud of cold, trapped 
atoms. 



2,1 Fundamental Approximations 

In the current generation of EEC experiments [1-3], 
a cloud of alkali atoms is optically pre-cooled and then 
magnetically trapped and evaporatively cooled to very 
low temperatures. The first major approximation leading 
to the MFT description is that the internal states of the 
atoms are ignored. All of the atoms must, however, 
reside in a particular hyperfine atomic ground state in 
order to remain trapped. The direction of the magnetic 
moment associated with the atom's internal state has 
been polarized to lie along the direction of the trap 
magnetic field at the site of the atom. Since the atoms 
are very cold and thus slowly moving, we assume that 
the magnetic moment of the atom adiabatically follows 
the local magnetic field [17]. Thus the energy of the 
interaction of the atom's magnetic moment (fi^tom) with 
the external magnetic field has the form 



length, and the interaction potential may be written in 
the form: 



Vtrap(/') = Atatom^(r). 



(1) 



Another feature of this assumption is that collisions be- 
tween atoms in the cloud do not change the atom's 
internal state. That is, all collisions are assumed to be 
elastic. In fact, most inelastic (spin-flip) binary colli- 
sions will cause both atoms to be ejected from the trap. 
This, in turn, limits the lifetime of the condensate. Such 
lifetimes can be predicted within MFT in a reasonably 
accurate way for comparison with experiment. Such 
comparisons are presented below. 

The true interaction potential between atoms in the 
cloud is quite complex. See, in this regard, Refs. [12-14] 
in this Special Issue. Most of this complexity is evident, 
however, only when the atoms are in close proximity. At 
the low temperature and density conditions present in 
the trap, all scattering events occur at extremely low 
energy. Consequently, the atoms rarely come close 
enough to each other to sample the complex nature of 
the inter-atomic potential. The atom-atom interaction is 
therefore well characterized by the ^-wave scattering 



VUr-n=UoHr-n 



(2) 



where Uq = Anfi^alM, a is the ^-wave triplet scattering 
length, and M is the atomic mass. 

In the next section we present the derivation of the 
MFT equations using the Eogoliubov prescription 
which begins with the assumption that the atomic cloud 
can be approximated by a restricted grand-canonical 
ensemble. 



2.2 Bogoliubov Prescription 

Consider the many-atom system whose temperature 
is well below the condensation point and which is com- 
posed of a condensate plus thermal atoms. The grand 
canonical, many-atom hamiltonian, K = H — [jlN where 
H is the many-body hamiltonian and N is the number 
operator, is written in terms of the field operator as 
follows: 

K=jdr ip\r)Hoip{r) + \ ^ojdr ip\r)ip\r)ip{r)ip{r) 
- fxjdr ip\r)ip{r) (3) 



where Hq is the bare-trap hamiltonian, 

Ho=-^V^+Vtrap(rX 



(4) 



IJL is the chemical potential, and VtrapC'') is the trap poten- 
tial. 

The boson field operators ip'^ir) and \p{r), respec- 
tively create and destroy an atom at position r and satisfy 
the commutation relations. 



[^(r),^V)] = W - r'\ 



(5) 



Under the Eogoliubov approximation, the condensate 
is assumed to contain most of the atoms so that 
N — No« No, where A^o denotes the macroscopic occu- 
pation of the condensate and N denotes the total number 
of condensate plus thermal atoms. In this case, the field 
operator can be written as the sum of a c -number con- 
densate wave function, '^(r), plus a small correction. 



iP{r)=nr) + 4>{r\ 



(6) 
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where ^(r) satisfies the normalization condition 

Idr l^(r)P=iVo. (7) 



h 



Inserting Eq. (6) into Eq. (3) and neglecting terms in 
(^(r) higher than quadratic yields the following expres- 
sion for K. 



K-- 



dr ^\r)\_Ho - At + 2 t/ol^(r)P]^(r) 



+ J dr ^»[//o - ju + t/oli/'(/-)P]^(r) 

+ J d/- ^\r)[H^ - (1+ [/ol^(r)P]^(r) 

+ J dr <pKr)[Ho - fj. + 2Uomr)\']4>(r) 

+ \ t/ojdr ^\r\^{r-)f^\r-) 

The first term in the above equation is a c -number and 
the second and third terms will vanish identically if 
^(r) satisfies the GP equation [18] 



[//o+ t/ol^(r)P]^(r) = ii^{r\ 



(8) 



The Bogoliubov-approximate grand canonical hamilto- 
nian [15], K^, then takes the form 

Zb = r + Jdr ^V)[^o - At + im^{rt\^{r) 
+ it/oJdrc^V)(^(r))W) 
+ it/oJdr^(r)(^V))^^(r), (9) 

where ^' is a c -number. 

The Bogoliubov hamiltonian is a sum of a quadratic 
form and a c -number and can be cast into the form of a 
collection of noninteracting quasi-particles by the fol- 
lowing Bogoliubov transformation [19] 

^(0 = S(^A(r)i8A + v;(Oi8l) (10) 



and 



^V) = S(^I(r)i8l + VA(r)i8A) (11) 



condensate wave function is not included in the sum. 
The quasi-particle operators satisfy the usual commuta- 
tion relations for boson creation and destruction opera- 
tors 

V^xm = V, Wx,M = WlPl] = 0. (12) 

The reduction of Kb to a collection of noninteracting 
quasi-particles occurs if the ux and vx satisfy the follow- 
ing equations (after setting ^(r) =No^ij/g(r)) 

Su,(r)-^NoUo(ips(r)fvx(r)=E,u,(r) (13) 

and 

Sv,(r) + NoUo(il^;(r)fu,(r) = - E,v,(r) (14) 

where 

2 = Ho- fJi-\-2NoUo\il^g(r)\\ (15) 

and the ux and vx are square-integrable functions. The 
final form (to within a c -number) of Kb is 



Kb = ^E,I31I3,. 



(16) 



where the Px are quasi-particle creation and destruction 
operators and the implicit assumption is made that the 



This hamiltonian has the form of a collection of nonin- 
teracting quasi-particles for which the condensate is the 
vacuum. Complete details of the derivation of the final 
form of Kb are given in Ref. [19] 

2.3 Linear-Response Theory 

The time-dependent Gross-Pitaevskii equation can be 
derived from an action principle if a ' 'boson coherent 
state" is used as the trial wave function [20], 

i«^(o)=i(§)"W„(o). (17) 

The factor l^„(0) is the usual Hartree many-body trial 
wave function — an n -fold product of one single-particle 
orbital ^(r,t). The time-dependent GP equation de- 
scribes the evolution of this orbital, 



i^^= [Ho^Uo\nr,t)P]nr,t). (18) 

at 



The basic MPT equations can be obtained by a standard 
linear-response analysis of this equation [16]. To this 
end, we consider the effect of adding a weak, sinusoidal 
perturbation to the trap potential. The time-dependent 
GP equation then has the form 
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ih^= \Ho+Uo\^{r,t)\^+U{r)c- 



dt 



+/-(r)e^"pT^(r,0. 



(19) 



The/+(r) are the (possibly spatially dependent) ampli- 
tudes of the sinusoidal perturbation and Wp is the probe 
frequency. 

To find the linear response of the condensate to the 
driving field, we shall assume that ^(r,t) takes the form 
of a sum of an undisturbed ground-state part and a 
response part that oscillates at frequencies ±0)^: 

^(r,t) = c-'^''\Nl ife(r) -F M(r)e-'"p^ -\- v(r)c''^'] (20) 

Here, fju is interpreted as the chemical potential of the 
undisturbed ground state, and the condensate wavefunc- 
tion is represented by the (scaled) condensate orbital 
ij/gir). The functions u(r) and v(r) are the components 
of the condensate's linear response to the external dis- 
turbance that oscillate at frequencies ±Wp. 

After inserting Eq. (20) into Eq. (19), retaining only 
terms up to first-order in u (r), v(r), and/± and equating 
like powers of e'^'^p^, there result three equations that 
must be simultaneously solved for ij/gir), u(r), v(r), and 
fjL. These equations describe the linear response of the 
condensate to weak external perturbation and have the 
following form 

[Ho^NoUo\il^Art]^s(r) = fJi^A^l (21) 

[S - ^co,Mr) + NoUo(il^Ar))Mr) = - Nl%(r)iljg(r\ 

(22) 

{2+fioyMr)^NoUo{ifj;{r))Mr) = - NiJ-(r)ilj;(r). 

(23) 

To make the final connection with the basic equations of 
MET, we show that Eqs. (22) and (23) can be solved by 
writing their solution as an expansion in the condensate 
normal modes. The equations that determine these 
modes are identical to the Bogoliubov equations. We 
now discuss there solution. 

To find the normal modes of the condensate, we first 
set/±(r) to zero in Eqs. (22) and (23). It is clear that the 
resulting equations will support square-integrable solu- 
tions only for discrete values of Wp, (we shall label them 
as (Ox). The normal- mode equations thus have the form 

[£ - ^co,]u,(r) + NoUo(iljAr)fyx(r) = 0, (24) 



and 



where A represents a set of quantum numbers. These 
equations are identical to Eqs. (13) and (14) if £'a = ^cox. 
To complete the connection between the quasi-particle 
excitation spectrum and the condensate response we 
now show how these normal modes describe the con- 
densate linear response. 

We define a normal mode as the following two-com- 
ponent object: 



Mr) 



/«A(r)\ 
'Wr))- 



(26) 



With this definition, Eqs. (24) and (25) can be cast in the 
form. 



H(l)x{r) = ^(Ox(T3(l)x{r). 



(27) 



Where 



with V(r) = NoUo(il/gf(r) and S is defined by Eq. (15). 
The {4>x} form a complete [20] orthonormal set 
where the scalar product of two normal modes is de- 
fined by 

<<^,J0,,) ^ Jdr 0l/r)cr30A,(r) = 8x,x, (29) 

and the t denotes the transposed, complex-conjugated 
matrix. 

The linear response equations, Eqs. (22) and (23), can 
be written, using this notation, as 



(H - ^(oa^jij/ir) = - a^girX 



(30) 



where 



NoUo(il^;(r)fux(r) + [S -^ ^cox]vx(r) = 0, (25) 



^(r)-(^^'A a(r)-( ^o%{r)il^g{r) \ 

The solution of the linear response equations is found by 
expanding both ijj(r) and g(r) in the normal modes 

il^{r) = 2cxMrX g{r) = ^gxcl>x{ry (32) 

A A 

Where the ^a are given by the following overlap integral 

gx = jdrcl>Kr)a,g(ry (33) 

Substituting these expansions into Eq. (30) yields a sys- 
tem of completely uncoupled equations to be solved for 
the cx . The final solution is written as 
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^(r) = - 2 



(gA^) 



&>A 



Mr)- 



(34) 



Note that the linear response diverges when the conden- 
sate is driven exactly on resonance. This unphysical 
behavior results from our neglect of loss processes and 
nonlinear effects. 



Mean-Field Theory: 
Experiment 



Comparison with 



Mean-field theory provides predictions for a variety 
of measurable condensate properties. These properties 
include condensate geometries, densities, and excitation 
frequencies. In this section, we shall discuss the com- 
parison of the predicted values of these quantities with 
experiment. Measurements of condensate properties 
have been, to date, primarily performed on condensates 
formed in the traps of Refs. [1] and [3]. Comparisons 
with properties of condensates formed in the trap of 
Ref. [2] are not possible at present because measure- 
ments of these properties still have substantial uncer- 
tainties. 

The geometry, density, and excitation frequency pre- 
dictions of the MFT equations do not contain any ad- 
justable parameters. The numerical constants that are 
input into the theory are the atomic mass (M), the radial 
and axial trap frequencies (ojp and w^.), the number of 
condensate atoms (A^o), and the scattering length (a). All 
of these numbers are determined experimentally. The 
scattering length, in particular, is determined princi- 
pally by photoassociation spectroscopic measurements 
[14,21]. 

Condensate lifetimes, on the other hand, depend crit- 
ically on two- and three-body scattering event rates. 
These rates are quite difficult to determine accurately 
[13]. Therefore we will not consider lifetimes in the 
comparison of MFT with experiment. 

3. 1. Geometries and Densities 

3.1.1 Geometries and densities were among the 
first condensate properties to be measured (see e.g., Ref. 
[1]). It is important to note that both of these ^^Rb 
condensate properties are determined indirectly. That is, 
after condensate formation, the trap potential is dramat- 
ically lowered, allowing the condensate to undergo a 
ballistic expansion, and then an absorption picture is 
taken. This picture, which essentially exhibits the veloc- 
ity distribution of the original condensate, can then be 
used to extrapolate back to the original spatial density 
distribution. 



The most quantitative comparison of theoretical and 
experimental geometry and density performed to date is 
that of Holland and Cooper (Ref. [22]). In this work, the 
time-dependent GP equation was solved by direct nu- 
merical integration. The time variations of the trap po- 
tential that occurred in the experiment of Ref. [1] were 
modeled and a prediction for the absorption picture was 
obtained. The agreement between theory and experi- 
ment for the velocity distribution appears to be at the 
5 % level. 

To date, no measurements of condensate geometries 
and densities of the ^^Rb condensate have been obtained 
in situ . However, it is interesting to see what the MFT 
predictions are for this case. Figure 1 shows the density 
profile for the ^^Rb condensate for the trap parameters of 
Ref. [1] for the case of A^o = 2012 atoms. A plot of peak 
density as a function of condensate population is shown 
in Fig. 2. The peak density for this case is 5. 2 X 10^^ 
atoms/cm^. The extrapolated experimental result is 
3 X 10'^ atoms /cm^ but this number is accompanied by 
substantial error bars within which the MVT number 
falls. 




z(|am) 



Fig. 1. A plot of the spatial distribution of the ground state ^^Rb 
condensate density over a plane containing the z-axis. In this plot, 
A^o = 2012 atoms. 




500 



1000 1500 
N (Atoms) 



Fig. 2. A plot of the peak condensate density as a function of conden- 
sate population for a ^^Rb condensate confined in the strong TOP trap. 
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3.1.2 ^^Na To date, condensates have been formed 
in two different traps at MIT. In the optical-plug trap, 
atoms were prevented from undergoing majorana spin 
flips through the use of a blue laser focussed at the 
center of the trap [3]. The "cloverleaf trap is an loffe- 
Pritchard trap and thus there is no zero of the magnetic 
field [23]. The condensates formed in each of these 
traps are very different from the ^^Rb one. The major 
difference is size — the ^^Na condensates contain on the 
order of 10^ atoms. 

Mean-field theory solutions for condensates of this 
size can be obtained via the "Thomas-Fermi" approxi- 
mation [24]. The approximation amounts to the neglect 
of the kinetic energy term in the Gross-Pitaevskii equa- 
tion. This reduces Eq. (8) to 



(Vtrap(r) +iVo^o| ^(r)p^(r) = fj.ijj(rX 



(35) 



where the kinetic energy term, - (^^/2M)V^i/^(r), has 
been neglected. The wavefunction thus has the form 



^(r)= i 



M-ytrap(^) 

NoUo 
0, 



for ytrap(^) < M ; 

fory,,,p(r)^)a, (36) 



The relationship between A^o and fjb is found by the 
normalization condition 



I 



dr\ilj(r)\^=\. 



(37) 



A comparison of the Thomas-Fermi approximate solu- 
tion with the basis-set solution is shown in Fig. 3. 
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Fig 3. Comparison of the Thomas-Fermi and MFT solutions for the 
MIT condensates. The condensate consists of 5 million ^^Na atoms 
confined in a cylindrically symmetric trap with Vp = 350 Hz and 
Vz= l^ Hz. The basis-set calculation required 350 functions. The plot 
shows the condensate density in the z = plane. Note that the correc- 
tion to the Thomas-Fermi result (shown inset) is important only at the 
very edge of the condensate. 



The comparison of the geometries and densities of 
the ^^Na condensates with the results of MFT suffers 
from the uncertainty in the value of the sodium scatter- 
ing length. In contrast to the ^^Rb scattering length, until 
recently the ^^Na scattering length was known only to 
within a factor of two [14]. The MFT solutions depend 
only on the parameter Noa/kho where /sho is the har- 
monic oscillator length scale. Uncertainty in a will lead 
to a similar uncertainty in MFT predictions. 

3.2 Excitations 

The excitation data provide the opportunity for the 
most quantitative comparison with MFT to date. Excita- 
tion frequencies are measured by forming a condensate, 
driving it weakly by oscillating the trap potential, wait- 
ing a specified delay time, and then probing the conden- 
sate. This cycle was repeated for increasing delay times 
forming a time history of condensate oscillations. The 
width of the condensate was observed to oscillate and 
the frequency of this oscillation measured. This experi- 
ment was recently performed both at JILA [25] and at 
MIT [26] . The comparison [27] with the JILA results is 
shown in Fig. 4. The agreement varies between 2 % and 
5 %. Similar agreement with the theory [28] was found 
in the MIT experiment. 



2.0 



I 1.8 



1.6 



1.4 



o 

B 1.2 



1.0 




Fig. 4. Comparison of the JILA excitations results with MFT predic- 
tions. This graph is reprinted from Edwards et al., Phys. Rev. Lett. 77, 
1671 (1996). 



3.3 Summary 

This paper has been written in the early days of quan- 
titative modelling of dilute atomic BECs, and there have 
yet been relatively few stringent tests of the validity of 
MFT. However, it has been found to have good predic- 
tive and interpretive value. As the accuracy of the exper- 
iments improves, and the uncertainties in the values of 
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the microscopic parameters are reduced, we expect that 
dilute atomic BECs will provide a new testing ground 
for the Bogoliubov approximation and its variants, 
which are among the cornerstones of the quantum the- 
ory of many-particle systems. 



4. Appendix A. Numerical Techniques 
4.1 Solving the Gross-Pitaevskii Equation 

The time- independent GP equation is written 

'2M 



V'+ytrap(r)+ NoUomr)V] ip(r) = fJiip(r). 



(38) 

We expand ipir) in the basis of trap-potential eigenfunc- 
tions, i.e.. 



ip(r) = 2 ax(i>x{r\ 



(39) 



where the sum extends over the A^basis basis-set functions, 
A represents a set of quantum numbers, and <l>\{r) is an 
eigenfunction of the trap Hamiltonian, 

HoMr) = {-|^ V^ + y,.p(r)} (/>,(r) (40) 



= £x(f>x{ry 



(41) 



Inserting Eq. (39) into Eq. (38) and taking the scalar 
product with each basis- set function converts the prob- 
lem to one of finding the simultaneous root of A^basis 
nonlinear algebraic functions having the following form 

f^{[a}) = (sa -l^)ak 
+ N0U0 2 C(A, Ai, A2, A3)aI^aA2«A3 = 0. (42) 

A J, A 2, A3 

The coefficients C(Ai, A2, A3, A4) are the integrals of the 
product of four trap wavefunctions; several methods for 
evaluating these integrals are described in Appendix 5. 
These nonlinear equations can be solved by a standard 
Newton-Raphson [29] technique provided a good guess 
at the solution is known. 

Equations (42) are not sufficient to determine the set 
{a } since ix is also an unknown. If we choose a value of 
IX , and a corresponding value of A^o, a solution may be 
found that will, in general, not be normalized. This 
solution may be normalized by the simultaneous rescal- 



ing of the basis-set coefficients {a}, and the number of 
condensate atoms, A^o- 

For example, find an arbitrary solution iI/na(^) of the 
time independent GP equation 



2M 



V^+ytrap(r)+A^. 



.Uo\lpNAr)\''j^N, 



(r) 



= fxij/MArX 

then find the norm of the solution 



I dr \il,^Ar)\\ 



(43) 



(44) 



where A is some constant. The set {a^}, the basis- set 
coefficients for the solution \fjM^{i'), and the condensate 
population Na may be rescaled as follows 



{«b} 



"UJ 



, andA^B = A'A^^, 



(45) 



to give a normalized solution i^Jn^T) for the new GP 
equation 

= lJi^N,ir). (46) 

So we find an arbitrary solution ij/na(^)^ then find a 
normalized solution ij/nb(^) from ij/na(^) t>y rescaling 
both the wavefunction and the number of atoms in the 
condensate. This method is also discussed in Refs. [30, 
31]. 

To ensure that the final normalized solution is actu- 
ally the lowest energy state of the condensate, we simu- 
late the adiabatic growth of the nonlinear term in Eq. 
(38) by first finding a solution for a value of fx close to 
the noninteracting ground state energy, which corre- 
sponds to a very small number of particles in the con- 
densate, and so the condensate wavefunction is very 
close to that of the trap ground state. The initial values 
of the basis-set coefficients are those of the trap ground 
state, i.e., {a} ={1,0,0, . . .}. The Newton-Raphson al- 
gorithm is then used to solve Eqs. (42). The chemical 
potential fx is then incremented, corresponding to a in- 
crease in condensate size, and the basis-set coefficients 
found previously, for the smaller condensate, are used as 
the initial values for this new solution. This process, 
which is repeated for increasing values of fx , gives not 
only the lowest energy solution, but also generates solu- 
tions for a large range of fx , and corresponding A^o, with 
little wasted computational time. 



560 



Volume 101, Number 4, July-August 1996 

Journal of Research of the National Institute of Standards and Technology 



4.2 Solving the Bogoliubov Equations 

In this section we describe how Eqs. (13) and (14) are 
numerically solved using a finite basis consisting of trap 
eigenfunctions. The solution is found by first expanding 
u\(r) and v\(r), in the eigenstates of Hq to produce a set 
of algebraic equations. These equations are then cast 
into a generalized matrix eigenvalue equation having the 
form Ava = w\Bx\. In this equation, A and B are square 
matrices of order 2A^basis (where A^basis denotes the num- 
ber of basis-set states used), X\ is an 2A^basis-component 
column vector containing the coefficients in the basis- 
set expansion of ux(r) and v\(r), and w\ is the general- 
ized eigenvalue. We shall treat the case of a spherically 
symmetric harmonic trap potential for simplicity. It is 
important when solving the normal-mode equations us- 
ing a finite basis that Eq. (21) for ijjgir) also be solved 
within this set. 

The exact expansion of ux(r) and v\(r) in eigenfunc- 
tions of Ho has the form 



Ux 



(r)= E 






^2, 



.ir) 



and 



where 



s\r,m' 



(47) 



(48) 



where /jl = jB^cotmp and cox = oix^trsL^- The full set of 
equations that must be solved to find the aA consists of 
those given just above for all values of / and m . The 
matrix element w^^?^ is given by 



w: 



ilm) , 



1 



^^rap 



f 



drw,/(r) 



NoUo\ip,{r)\^ 



Usi(r). (52) 



where \p2s^i,i,mi^) = Usi(r)Yimif)lr defines the reduced 
radial wavefunction Usi(r). 

The spherical symmetry of the condensate ground 
state imposes a particularly simple structure upon the 
equations for ashl and Z?i£. These equations divide into 
subsets of equations that are self-contained within the 
subspace of fixed angular momentum and magnetic 
quantum number, i.e., there is no coupling between a^l'^ 
or b^tJi and either aifL' or bifl^'. Thus an infinite number 
of normal-mode solutions, each proportional to a single 
Yim , can be found by setting to zero the a and b coeffi- 
cients corresponding to all pairs (/', m') except one [say 
(/, m)] in expansion Eqs. (47) and (48). We can replace 
the eigenvalue label A with the quantum numbers / and 
m plus a "radial" quantum number n. The significance 
of this quantum number can be understood by writing 
the normal modes in terms of their basis-set expansions 
with the above symmetries imposed 



^basis-l 



Unim{r) = -YUr) 2 a^lu,i{r)^UnUr)YUr){S3) 



Hollj2s + U,m(r) = (2s -F Z-F 3/2)MraptA2. + /,/,m(r) 

^4^->MraptA2../,/,.(r) (49) 

defines the eigenvalues and eigenfunctions of //sho- The 
ranges of variation of the quantum numbers are ^ , / = 0, 
1, . . . and m = -/,..., /. 

Truncating and inserting these expansions into Eqs. 
(24) and (25) and projecting each equation onto an arbi- 
trary eigenstate of Hq yields the following system of 
linear equations 

A^basis-l 



+ E w^?"^it = (50) 

^• = 



and 



and 






(A) 



A'tosis-l 

+ 2 2 yv'i?blt = 0. (51) 



^basis-1 



Vn,n,{r) = - Y,^(f) Si b^^ uAr) ^ V„,^(r)Y,^(r). (54) 

^ s = 

The quantum number n is defined as the number of 
nodes found in the radial wavefunctions, Unim{r) and 
Vnimir). Furthermore, n runs from -oo to co because the 
normal frequencies occur in pairs of opposite sign. That 
this must be the case can be understood by examining 
Eqs. (24) and (25). If a>A^ - ojx and if ux and vx are 
exchanged, then Eqs. (24) and (25) remain unchanged. 
Hence, if there exists a solution of the normal-mode 
equations for positive wa, then a negative frequency 
solution also exists having the roles of ux and vx reversed. 
Indeed, we find that normal frequencies occur in equal- 
magnitude pairs of opposite sign in our solutions. A 
negative value of n, therefore, labels a normal-mode 
solution obtained by performing the above transforma- 
tion on the corresponding positive-n solution. The n = I 
= m = normal-mode solution deserves special com- 
ment. This solution always has a frequency of zero (i.e., 
a>ooo = 0) with Wooo('') and Vooo('') being proportional to 
t/fg(r). This condensate-mode solution has zero norm, 
see Eq. (29). 
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Moving the terms in Eqs. (50) and (51) containing a^ 
to the right-hand- side gives 



A^basis-l 



A^basis-1 



(^r^-js)ai£ + 2 2 wr^i?i+ 2^^^^. 






(55) 



and 



(^<"">-;S)Z.«+ 5 wr«i"i + 2 E w^^i"; 






(56) 



where we have replaced ax with a^M and (a^ti, Z^it^i ) 
with aiti, b^^li ). The above A^basis pairs of equations can 
be summarized as a single (2 X 2) matrix equation by 
making the following definitions. 



, nlm _ I ^slm 



n>h (^■■ 



2 1 



, cr3 



1 



(57) 



Kb^S./ ' " \1 2/ ' "^ VO -1, 

With these definitions, Eqs. (55) and (56) can be written 

as 

s' = 

(58) 

where /s is the (2 X 2) unit matrix and where s runs 
from to A^basis-1- 

The final transition to the form of a generalized ei- 
genvalue problem is accomplished by defining the fol- 
lowing (A^basis X A^basis) matHccs 



and 



(d'^-\, = (^i^-^ - 13)8,,, 



(w^'-'X' = wi-\ 



(59) 



(60) 



Using these matrices, Eq. (58) may be written as a 
single matrix equation 



j^(lm)^(nlm) _ ^ ^(Im) ^(nlm) 



(61) 



where c ^" ^^ is given by 



, (nlm) _ 



, (nlm) 



'^A^ basis- 1 



and where the matrix A ^^""^ can be expressed as a tensor 
product of matrices: 

and B is given by 

B""' = I^^ ® ^3. (64) 

where iN^^^i, is the unit matrix of order A^basis- 

The square matrix on the left of the ® symbol is of 
order A^basis and the matrix on its right is of order 2. The 
full matrix is of order 2A^basis with elements that can be 
indexed as C(sMs'k')='^ss'Bkk', where s runs from to 
A^basis - 1 and k runs over 0, 1 . 

The component expressions for the matrices appear- 
ing in the final generalized-eigenvalue problem are as 
follows 



(A'^isMsT) = (d'nss'8kk'-^(w''nsAl + 5,,0 (65) 



and 



(66) 



(^ ^ ) (sMs'k') — Oss'Okk' (OkO— Oki). 

Furthermore, 

(nlm) _ ^ (n) « (nlm) _ i (n) ff^l\ 

C(s,0) —CI slm, ^(s,V) — ^ slm- V^ ' ) 

5. Appendix B. Integration Methods 

The basis-set expansion method requires the evalua- 
tion of a set of integrals, the integrals of four basis- set 
functions over all space. In general the basis-set func- 
tions would be those for an anisotropic trapping poten- 
tial, though in two of the cases of interest [1,2] they are 
eigenfunctions of the cylindrically symmetric harmonic 
oscillator. We will now consider how these integrals 
may be evaluated. 

5.1 Cylindrical Symmetric Traps 

The integral which must be evaluated is given by 

C(Ai, A2, A3, A4)= d(^ dpp 

I Az(i>\Ar)(i>\^{r)(i>,,{r)(i>,,{r) (68) 

where the basis-set functions <f>\{r), for a harmonic 
trapping potential 



(62) 
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are given by 



(69) 



0A,(r) = 



175 2 "^'T?,.! 



e-5"'^^'H„.(a,z) 



X 



'n(np. + mi)\ 



(70) 



where p = (x^ + yy^ and dp,, = (Mwp,,/^)"l These 
functions may be rewritten, with substitutions 
^= (appf and ^= a,z, to give 

ct>,Xr) = N,, e-5f'H„^_. (O e""-* f? e^^ L<„-;;>(D, (VI) 
the normalization constant remains unchanged, 



A^A 



^pfip.l 



and the integral operator transforms, 



(72) 



(73) 
These substitutions reduce Eq. (68) to 



C(Ai,A2,A3,A4) = ■ 



M^3A^r 



d<h Q^^(^^+^4-fni-m2) 



•^0 



xLi.7;(OLi.7;(D 

X f d^ e-^f^ H„,, (^) H„,^ (^) H„,3 (^) H„,^ (^), 

(74) 

so the integral has been broken down into three one- 
dimensional integrals, i.e., 

C(Ai,A2,A3,A4) = — 2 — ' ^ (Ai, A2, A3, A4) 

2aoa, 



Xi7(Ai, A2, A3, A4) H (/2zj, n,^, ?i,3, n,^) (75) 



where 



0(Ai,A2, A3, A4): 



^{m^+m^-mY-m'i) 



(76) 



dc^ e''^^" 

U{\u A2, A3, A4) = d^ ^-2^ ^i(mi+;n2+-3+-4) 

X L L?V^) /^ L':?(^) L <7Y^) L ^-;KD (77) 



"P2^^^ "P3 



and 



f d^ e-^^X/^)^«.2(^)^«.3 (^)^«.4 (^)- (78) 

J —00 

5.2 Completely Anisotropic Traps 

The integral to be evaluated for a completely an- 
isotropic trap, 

Vtrap(r) = M(oj^x^ + oj^y^ -^oj} z^l, (79) 

is given by 

A(Ai, A2, A3, A4)= dx dj 

J -00 J —(X> 

/»C0 

J dz<t>\,{r)4>\,{r)<t>,,{r)4>,,{r) (80) 

where the basis-set functions c^a('') are now given by: 



<i>^i{r) = 



7ri2"-'7i,.! 



X 



X 



7ri2">'7i^.! 



a. 



7ri2"^'7i,.! 



^--"^"H^^ia^x) 



c-l-yy' H^^^(ayy) 



c-\^"'" H.{a,z) 



(81) 



where q;;,,^,,^ = (Mw;,,3.,2/^)^^^These functions may be 
written. 



X H^^,(a^x) Hn^,(ayy) H^^,(a,z). 



(82) 



With the substitutions axX = u,ayy = v and azZ = w, 
this reduces Eq. (80) to 
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A(Ai, A2, A3, A4): 



Nx,NxJ^x,Nx, 



a^^aya^ 



Xj dwe-^"'//„^^(w) //„,^(w) //„,3(w) H,,^(u) 
X J du e-^-'//„^^ (v) H,^^ {v) H^^^ {v) H^^^ {v) 

X J dwe-^-'//„.^(w) H,^^{w) H^^^iw) Hr^^^{w) 

(83) 



A(Ai, A2, A3, A4) = 



Nx,Nx,Nx,Nx, 



X B(n^^, n^^, n^^, n^^Biriy^, riy^, Uy^, Uy^ 

X S(n,^,n,^,n,^,n,J (84) 

we will now outline some of the methods we may use to 
evaluate these integrals. 

5.3 Exact and Series Methods 

Equation (76) may be evaluated easily to give 



277, for mi+m2=m3+m4; 



^(Ai,A2, A3, A4) = - 



0, otherwise . 



(85) 



Equation (78) may be evaluated by use of an addition 
formula for Hermite polynomials [32] to give 



'n3\n4l ™n^3,«4) 2^+^-1 r(s-ni) r(s-n2) 
^ h t\(n,-t)\(n,-ty. 

X r (s-n3-n4+2t), 

if ni+n2+n3+n4 is even; 



H(ni,n2,«3,«4)= ■< 



^0, 



otherwise. 



(86) 



where s = {n\+n2+n3+n4-2t+\)l2. 

By expanding the Laguerre polynomials in Eq. (77) 
as finite series [33], an even more complicated expres- 
sion can be derived for the cylindrically symmetric case. 
While these results are exact, they are of little practical 
use for computational calculations in standard extended 
precision arithmetic, since numerical accuracy is 
rapidly lost in the summations of terms which vary over 
several orders of magnitude. An alternative method for 



evaluating these integrals is that of gaussian quadratures 
which is described next. 

5.4 Gaussian Quadrature Methods 

For simplicity 

J d^e-^^'//„,(^) H^,{^) H^,(0 H^4(a 

(87) 

may be rewritten, with the transformation v= v2^, as 

H(/2i,/22, ^3, ^4) = 

Using a gaussian quadrature integration technique [29] 
reduces this integral to a finite sum 

H(/2i,/22, n3,n4) = 

1 ^Sad y (h) . ^ (h) . 

where g[^^ is the iih Gauss-Hermite quadrature point, 
and Wi^^ is the corresponding weight. This sum is 
exactly equal to the integral provided A^ [|]ad > 2 
max(?zj [29]. 

Applying the same technique to 

t/(Ai, A2, A3, A4) = f d^e-^^^I^-^^-^^-^^-^^ 

X L ':i\o L ':i\c) L ':i\c) L ':f{c) (90) 

with transformation w = 2^, gives the sum 

1 



t/(Ai, A2, A3, A4) - 



2 2 (mi-\-m2-\-m3-\-m4-\-3) 



X 



Quad / \ 1 / ^ (.'A 
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where gp is the ith Gauss-Laguerre quadrature point, 
and wp is the corresponding weight. This sum is exact 
provided NqLd > 2 max(/2p)+max(m) [29]. 

This means that the choice of basis-set functions fixes 
the number of quadrature points needed, e.g., for initial 
calculations of a cylindrically symmetric ground state 
condensate, we may use /2p = {0, . . . , 6}, m = 0, and 
/2, = {0, . . . , 6}, fixing TV ^^ad = 13, and N'i^Ld = 13. 

Acknowledgments 

This work was partly supported by the U.K. Engi- 
neering and Physical Sciences Research Council, and by 
the National Science Foundation grants PHY-9601261 
and PHY-9612728. 



6. References 



[1] 



[2] 



[3] 



[4] 
[5] 

[6] 



[7] 



[8] 



[9] 



[10] 
[11] 
[12] 

[13] 

[14] 



[15] 
[16] 

[17] 



[18] 



M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, 
and E. A. Cornell, Science 269, 198 (1995). 

C. C. Bradley, C. A. Saclcett, J. J. Tollett, and R. G. Hulet, Phys. 
Rev. Lett. 75, 1687 (1995). 

K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, 

D. S. Durfee, D. M. Kurn, and W. Ketterle, Phys. Rev. Lett. 75, 
3969 (1995). 

S. N. Bose, Z. Phys. 26, 178 (1924). 

A. Einstein, Sitzber. Kgl. Preuss. Akad. Wiss. XVIII-XXV, 3 

(1925). 

A. Griffin, Excitations in a Bose-Condensed Liquid, 

Cambridge University Press, Cambridge (1993); H. Glyde, 

Excitations in Liquid and Solid Helium, Oxford University 

Press, Oxford (1994). 

A. Griffin and S. Stringari, Phys. Rev. Lett. 76, 259 (1996). The 

dilute limit of "^He can also be attained in a porous medium. See 

J. D. Reppy, Physica 126B, 335 (1984); P B. Weichman, M. 

Rasolt, M. E. Fisher, and M. J. Stephen, Phys. Rev. B 33, 4632 

(1986). 

P A. Ruprecht, M. J. Holland, K. Burnett, and M. Edwards, 

Phys. Rev. A 51, 4704 (1995). 

R. J. Dodd, M. Edwards, C. J. WiUiams, C. W. Clark, M. 

Holland, P A. Ruprecht, and K. Burnett, Phys. Rev. A 54, 661 

(1996). 

H. Stoof, Phys. Rev. A 49, 3824 (1994). 

E. Cornell, J. Res. Natl. Inst. Stand. Technol. 101, 419 (1996). 
P S. Julienne, J. Res. Natl. Inst. Stand. Technol. 101, 487 
(1996). 

F. H. Mies, C. J. Williams, P. S. Julienne, and M. Krauss, 
J. Res. Natl. Inst. Stand. Technol. 101, 521 (1996). 

E. Tiesinga, C. J. Williams, P S. Julienne, K. M. Jones, P D. 
Lett, and W D. Phillips, J. Res. Natl. Inst. Stand. Technol. 101, 
505 (1996). 

N. Bogoliubov, J. Phys. (U.S.S.R.) 11, 23 (1947). 
P.-G. de Gennes, Superconductivity of Metals and Alloys, Ben- 
jamin, New York (1966). 

A variant of MFT can be formulated without the adiabatic 
assumption: see Tin-Lun Ho and V. B. Shenoy, Phys Rev. Lett, 
(to be published). 

V. L. Ginzburg and L. P Pitaevskii, Zh. Eksp. Teor. Fiz. 34, 
1240 (1958) [Sov. Phys. JETP 7, 858 (1958)]; E. P Gross, J. 
Math. Phys. 4, 195 (1963). 



[19] 



[20] 



[21] 



[22] 
[23] 

[24] 
[25] 

[26] 



[27] 

[28] 
[29] 



[30] 
[31] 

[32] 
[33] 



A. L. Fetter, Ann. of Phys. 70, 67 (1972); A. L. Fetter, Phys. 

Rev. A 53, 4245 (1996). 

J. -P. Blaizot and G. Ripka, Quantum Theory of Finite Systems, 

MIT Press, Cambridge, MA (1986). 

J. R. Gardner, R. A. Cline, J. D. Miller, D. J. Heinzen, H. M. 

J. M. Boesten, and B. J. Verhaar, Phys. Rev. Lett. 74, 3764 

(1995). 

M. J. Holland and J. Cooper, Phys. Rev. A 53, R1954 (1996). 

M.-O. Mewes, M. R. Andrews, N. J. vanDruten, D. M. Kurn, D. 

S. Durfee, and W Ketterle, Phys. Rev. Lett. 77, 416 (1996). 

G. Baym and C. Pethick, Phys. Rev. Lett. 76, 6 (1996). 

D. S. Jin, J. R. Ensher, M. R. Matthews, C. E. Wieman, and 

E. A. Cornell, Phys. Rev. Lett. 77, 420 (1996). 

M.-O. Mewes, M. R. Andrews, D. M. Kurn, D. S. Durfee, 
C. G. Townsend, and W Ketterle, Phys. Rev. Lett. 77, 988 
(1996). 

M. Edwards, P A. Ruprecht, K. Burnett, R. J. Dodd, and C. W 
Clark, Phys. Rev. Lett. 77, 1671 (1996). 
S. Stringari, Phys. Rev. Lett. 77, 2360 (1996). 
W H. Press, S. A. Teukolsky, W T. Vetterling, and B. P 
Flannery, Numerical Recipes in C, Cambridge University 
Press, Cambridge (1992). 

Mark Edwards and K. Burnett, Phys. Rev. A 51, 1382 (1995). 
P. A. Ruprecht, K. Burnett, Mark Edwards, and C. W Clark, 
Phys. Rev. A (in press). 

L W Busbridge, J. London Math. Soc. 23, 135 (1948). 
M. Abramowitz and I. A. Stegun, eds.. Handbook of Mathemat- 
ical Functions, Dover, New York (1965). 



About the authors: Mark Edwards is an Associate 
Professor in the Department of Physics at Georgia 
Southern University. R. J. Dodd is a graduate student in 
the Chemical Physics program of the University of 
Maryland, College Park. Charles W. Clark is Chief of 
the Electron and Optical Physics Division of the NIST 
Physics Laboratory. K. Burnett is a Professor of Physics 
at the Clarendon Labratory of Oxford University, and is 
a Fellow of St. John's College. M. Edwards, R. J. Dodd, 
and K. Burnett are also guest researchers in the NIST 
Physics Laboratory. The National Institute of Standards 
and Technology is an agency of the Technology Adminis- 
tration, U.S. Department of Commerce. 



565 



